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ABSTRACT 


The objective of this talk is to describe three research thrusts in crashworthiness analysis: 

1) adaptivity 

2) mixed time integration, or subcycling, in which different timesteps are used for different parts of 
the mesh in explicit methods 

----- 3) methods for contact-impact which are highly vectorizable. 

The techniques are being developed to improve the accuracy of calculations, ease-of-use of 
crashworthiness programs and the speed of calculations. The latter is still of importance because 
crashworthiness calculations are often made with models of 20,000 to 50,000 elements using explicit time 
integration and require on the order of 20 to 100 hours on current supercomputers. 

The methodologies will be briefly reviewed and then some example calculations employing these 
methods will be described. The methods are also of value to other nonlinear transient computations. 
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OUTLINE 


• Adaptive mesh procedures in nonlinear 
analysis: why, how, and what is the 
status 

• Subcycling (mixed time integration) 

• New highly vectorizable methods for 
contact impact which are well suited to 
adaptive methods 

Figure 1 


PREDICTION 

The 1990’s will be the decade of adaptivity, 
adaptive mesh refinement 
adaptive targeting 
adaptive organization objectives 


Figure 2 



PREDICTION 


There are three types of adaptivity, which are known by the letters r, h, and p. These letters are 
mnemonic letters and refer to how the refinement is achieved. In r methods the nodes are relocated. In h 
methods, refinement is achieved by reducing the element size h. In p methods, refinement is achieved by 
increasing the order p of the element interpolance. 


TYPE OF MESH ADAPTIVITY 


r — method 



relocate nodes 


h — method 



element size h 


p — method 



adapt order p 


of element 


interpolants 


Figure 3 



ADAPTIVITY IN NONLINEAR FEM 


Adaptive methods are particularly useful in nonlinear problems such as crashworthiness because 
nonlinear response is often characterized by localization. In the areas of localized response more 
refinement is needed. When standard method is used, the user of the program must refine the mesh where 
he anticipates this localized deformation. Therefore, different meshes must be developed for different 
loadings. For example, in car crash, different meshes must be developed for frontal and rear impact, side 
impact, and overturning. This can be quite expensive from the viewpoint of manpower. 


Why are adaptive methods particularly important in nonlinear 
problems? 

Modes of failure of structures 

i. buckling, particularly with formation of hingelines 

ii. localization 

iii. fracture 


All of these involve local phenomena whose location cannot be 
determined at the outset of a simulation. 


Figure 4 



COMMENTS ON ADAPTIVITY FOR SHELL AND 
CRASHWORTHINESS PROBLEMS 


In comparing the different types of adaptivity for nonlinear structural dynamics problems such as 
crashworthiness, the following advantages, which are marked by a plus sign (+), and disadvantages 
which are marked by a minus sign (-), can be attributed to the various types of methods. From this study 
we concluded that the h-method was the most suitable method for adaptivity in crashworthiness. 


r — method 

- large elements cannot represent shape of shell 

+ most accuracy with given NDOF compared to h 

- history diffusion 

- elements become distorted - decr ease s accuracy 
+ easiest data structure 

p — method 

- awkward in nonlinear dynamics; no good lumped mass 
+ easy data structure 

h — method 

+ relatively effective 

+ no distortio n of elemen ts 

- moderately complex data structure 


Figure 5 
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TYPES OF ERROR INDICATORS 


Error indicators are an important ingredient in adaptive methods since they are to a guide the refinement 
of the mesh. Error indicators are classified by Oden in the following classes: residual, interpolation, and 
post-processing. In the work we are doing, we are using projection error criterions, a post-processing 
type, because they are very easy to implement and are quite effective for low-order elements. 


1. Residual: Compute residual in governing equations and use its 
norm or use it to drive an element or local enriched solution. 

a) Explicit: Evaluate a norm of the residual. 

b) Implicit: Use residual to drive a local or element error 
equation. 


2. Interpolation Methods: Estimate magnitude of derivatives of 
higher-order than contained in finite element space. 


3. Projection (postprocessing) Methods: Obtain a smoothed solution 
and compare to finite element solution; sometimes called L2 
projection methods. 


Figure 6 


14 



ADAPTIVE SCHEMES FOR TRANSIENT AND 
NONLINEAR PROBLEMS 


based on constant resource approach 

1 . Advance the solution n time steps 

2. Compute element error indicators 0 e 

3. Sort 0e 

4. Fission elements with 0 e > mission 
Fuse elements with 0 e < tol fusi ° n 

5. Repeat the last n time steps with new mesh (optional) 

6. go to 1 


Note: If n is too small or tol fission too close to tol fusion , we 
encounter "churning" which degraded accuracy. Our recent 
experience shows 5 is quite important. 


Figure 7 



REMARKS ON H-ADAPTIVITY 


Constraints (or slave nodes in explicit methods) must be introduced 
at nodes where a large element has two or more neighbors on one 
side to enforce compatibility; easy in vector methods, awkward in 
matrix methods. 


Usually a group of contiguous elements should be fissioned 
simultaneously because fissioning a single element does not provide 
much enrichment; only one new free node. 


In wave propagation problems, change in element size can cause 
spurious reflections. 


Usually mesh gradation is limited to 1 -irregular meshes: large 
element cannot have more than 2 small neighbors on any side; see 
Devloo, Oden and Strouboulis (1987). 


Data structure with fission and fusion is complex, particularly for 
real engineering meshes; see Belytschko, Wong and Plaskacz, 
Computers and Structures, 33(4-5), 1989, pp. 1307-1323. 


Figure 8 
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MIXED TIME INTEGRATION 


In h-adaptive meshes, large variety of element sizes are found. When explicit methods are used of 
such meshes, the timestep is reduced dramatically by the presence of small elements. Therefore methods 
called mixed time integration (or subcycling) are being u sed. 


Motivation : in explicit integration with same At over 
entire mesh, stiffest element sets At. also called subcycling, 
explicit-explicit partitions; 


example 










10 — *1 1 hb- 

LN 


— H 

h H— 

/ 


B A 


Atcrit = min (“) 
for A 

for AuB 


c = wave speed 



AWit- iQ c 


so AuB is lOx as expensive as A 


Figure 9 
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Mixed Time Integration 


Integrate each element or subdomain with At cr it using an interface 
treatment that preserves stability + consistency. 


In example 


2 4 



1 3 


10 — H I h h — 


integrate element 1 and nodes 1 to 4 with 



elements 2 to 10 and remaining nodes with 
At.i 

cost savings: ~ 90% 

In adaptive methods, large range of stable time steps is unavoidable, 
so subcycling is crucial for efficiency. 


Figure 10 
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CONTACT-IMPACT 


The modeling of contact-impact is very important in the simulation of crashworthiness. However, 
contact-impact algorithms often require more than fifty percent of the running time of a crashworthiness 
code because they are not easily vectorized. Therefore we have developed a pinball algorithm which is far 
more highly vectorizable. 


Contact-impact is an important phenomenon in crash analysis, e.g., 

1. engine impact with body, fire wall 

2. wheel impact with inner fender 

3. contact of collapsing surfaces 

Most contact- impact algorithms require many different branches. 



The branch of the algorithm which is activitated depends on which 
surface is penetrated; there are special branches for edges, etc. 


Figure 1 1 
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PINBALL PENALTY ALGORITHM 


T. Belytschko and M. O. Neal, International Journal for Numerical 
Methods in Engineering , 31, 1991, pp. 547-572. 


Interpenetration and inteij>enetration rate g are 
computed on pinballs inserted in elements. 



Enforces contact-impact conditions on spheres embedded 
in elements. 

As h — > 0, impenetrability is enforced. 

Algorithm is simple and highly vectorizable. 


Figure 12 
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Salient Features of Algorithm 


Radius of pinball is determined by equivoluminal expression 




47C 


Pinballs are classed by body; for single-surface slideline, smaller R 
needed. 


Interpenetration has occurred when 

d. < R + R 
y 1 j 

g = d- 
& y 


Pinball forces are equally transferred to all nodes of associated 
element (a surface node option available). 


The pinball method automatically places pinballs on outside 
elements by using assembled surface normal algorithm. 


Figure 13 
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EXAMPLES OF NONLINEAR ADAPTIVE 
COMPUTATIONS 


Nonlinear, transient computations with an explicit nonlinear finite 
element program WHAMS using h-adaptivity and pinball for contact 
impact; see Belytschko and Yeh (1992). 


An L2 projection on the strain invariants was used to calculate an 
error estimate. 


A commercial version of this program is available from: 

KBS2, Inc. 

455 Frontage Road 
Burr Ridge, IL 60521 
(708) 850-9444 
Fax (708) 850-9455 


Figure 14 
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TWO-LEVEL ADAPTIVE MESH OF CYLINDRICAL PANEL 


This shows an h-adaptive solution of a cylindrical panel which is impulsively loaded. Notice that the 
elements are refined along the side and at the support, where there is severe plastic bending deformation, 
and hinge lines consequently form. 



Two-level adaptive mesh of cylindrical panel. 


Figure 15 
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This figure shows a comparison between solutions obtained by h-adaptivity and those obtained using a 
very fine mesh and a coarse mesh. As can be seen, the adaptive solution compares well to the fine mesh 
solution. The differences in the displacements obtained by the coarse mesh and the fine mesh are not 
large, but for some of the stresses and strains, significant improvement is obtained by the use of 
adaptivity. 



Tune 


Displacement of point A 
in cylindrical panel ( MAXLEV=2) 



Tone 



Time 


Strain Exx of element B 
in cylindrical panel (MAXLEV=2) 



Strain Eyy of element B Stress Sxx of element B 

in cylindrical panel (MAXLEV=2). in cylindrical panel (MAXLEV=2) 



Stress Syy of element B in cylindrical panel (MAXLEV=2). 


Figure 16 
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This shows results for an S-beam which is impulsively loaded. Again, significant differences occur in 
some of the strains for a coarse mesh solution as compared to an adaptive or fine mesh solution. 





2 " ^ 


n 

2 " 


0.0472 


H 


3.6" 


T - shape cross section 

Material number 2 ( Table 6 ) 

Geometry of T-shape cross section beam. 



Axial displacement Dz of point A 
in T-beam (MAXLEV=1). 



Strain Exx of element B 
in T-beam (MAXLEV=1). 



Axial velocity Vz of point A 
in T-beam (MAXLEV=1) 



Strain Eyy of element 
in T-beam (MAXLEV=1) 


Figure 17 
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This shows the evolution of the mesh for the S-beam; note that the h-refinement occurs at a comer 
where local buckling takes place. 





Time = 1.28 ms ; 673 elem 


Time=1.92 ms ; 688 clem 


One-level adaptive mesh of T-beam. 
Figure 18 
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This again compares displacement in strains for coarse mesh, fine mesh, and adaptive solutions. 
Again the adaptive solutions agree very closely with the fine mesh solution. 



Axial displacement Dx of point A 
in T-beam (MAXLEV=2). 



Strain Exx of element B 
in T-beam (MAXLEV=2) 



Stress Sxx of element B 
in T-beam (MAXLEV=2). 



Axial velocity Vz of point A 
in T-beam (MAXLEV=2) 



Strain Eyy of element B 
in T-beam (MAXLEV=2) 



Stress Syy of element B 
in T-beam (MAXLEV=2). 


Figure 20 





This is a solution of a box beam which has an initial velocity as shown with an attached mass at the 
back. This problem is often considered a model for crash analysis. The solutions for fine mesh, coarse 
mesh and adaptive meshes are shown; the adaptive solution agrees well with the fine mesh. 



Geometry: L = 0.1500 m Initial condition: V=15.64 m/sec 

a - 0.0300 m Material number 3 (Table 6) 

t = 0.0015 m 

Box beam problem. 



T«* 

Velocity of point A in the box-beam (MAXLEV=1). 



Reaction force of box-beam (MAXLEV=1). 
Figure 21 
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This shows a timing for a full car model which is shown on the next page. It is solved with full 
contact-impact and subcycling. The important thing to notice is that subcycling gives a speedup of 1 7 and 
that the effective element cycle time on a CRA Y- YMP here is 1 2 microseconds. 


Timing 

FULL-CAR MODEL 


Elements: 

17,297 

Mass (kg): 

1,880 

Time steps: 

78,274 

80 msec simulation 


CRAY- YMP 


Without subcycling: 
128 elements/block 

7.63 hrs 

Element cycle time: 

20 fisec 

With subcycling: 
64 elements/block 

4.39 hrs 

Effective element cycle time: 

1 2 jisec 

Speedup: 

1.7 
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Figure 22 



wlnBB - von-mises shell material - subcycle 
time ■ 0.000E+00 



Figure 23 
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uiinBQ — von-mises shell material - subcycle 
time ■ 0.000E+00 



Figure 24 


32 


)nu»itn%\v 


uiinBB 

time 




von-mises shell material - subcycle 

A . 001E+01 



Figure 25 
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Implementation and Stability 

Time steps are assigned to nodes and element blocks automatically. 


Elements are sorted by At£ nt 


At crit < At crit < 

e 2 


< At crit 

n 


Elements are arranged in blocks so that time steps of adjacent 
elements have integer ratios.* 


Blocking of elements is necessary to take advantage of vectorization. 

For analysis of stability, see Belytschko and Lu, ASME publication 
edited by G. Hulbert, et al, 1992. 


*A new algorithm which does not require integer ratios has recently 
been developed (Belytschko and Neal, Computer Methods in 
Applied Mechanics and Engineering , 31, 1989, pp. 547-570). 


Figure 26 
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REMARKS AND CONCLUSIONS 


• H-adaptivity is a promising technique for 
simulating nonlinear structural response and 
structural failure. 

• Improves accuracy. 

• Simplifies model preparation. 


• Subcycling and advanced contact-impact methods 
such as the pinball method can improve efficiency 
of explicit dynamic codes and is essential with h- 
adaptivity. 


• Improved error criteria are needed for adaptive 
methods for nonlinear solid mechanics. 


Figure 27 
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